Analysis of GPP response to N and P enrichment

Libraries

library(tidyverse)
library(gamm4)
library(mgcViz)
library(MuMIn)
library(forecast)
library(gridExtra)
library(tidymv)

Load & munge data

Loads data, calculates log GPP and ER, does a ridiculous transformation of NEP, and centers stream water temp at 10C as inverse kT

met2 <- read.csv("00_metab_for_Hoodie.csv") %>% 
  mutate(Pdt = as.POSIXct(Pdt, format = "%m/%d/%y"),
         doy = as.numeric(strftime(Pdt, format = "%j")),
         Y = as.numeric(strftime(Pdt, format = "%Y")),
         Yf = as.factor(as.character(Y)),
         lGPP = log(GPP_median),
         lER = log(-ER_median),
         #centering
         invKT.C = one_over_kT - 1/((10 + 273.15)*0.0000862),
         streamTreat = as.factor(paste0(stream,"_",treatment))) %>% 
  mutate_at(vars(stream, treatment), funs(as.factor(.))) %>% 
  rename(invKT = one_over_kT)

For predictor variables, we want to center Temp, light, Q, and GPP. We’ll do that using stream means, which we calculate below

# stream mean
metS <- met2 %>% 
  group_by(stream, treatment) %>% 
  summarise_at(vars(meanTemp, invKT.C, meanQds, LightPerDay, lGPP), funs(mean(., na.rm = TRUE))) %>% 
  group_by(stream) %>% 
  summarise_at(vars(meanTemp, invKT.C, meanQds, LightPerDay, lGPP), funs(mean(., na.rm = TRUE))) %>% 
  rename_at(vars(-stream), function(x) paste0(x,".StMean"))

# global mean
metG0 <- met2 %>% 
  group_by(stream, treatment) %>% 
  summarise_at(vars(meanTemp, invKT.C, meanQds, LightPerDay, lGPP), funs(mean(., na.rm = TRUE))) %>% 
  group_by(stream) %>% 
  summarise_at(vars(meanTemp, invKT.C, meanQds, LightPerDay, lGPP), funs(mean(., na.rm = TRUE))) %>% 
  group_by() %>% 
  summarise_at(vars(meanTemp, invKT.C, meanQds, LightPerDay, lGPP), funs(mean(., na.rm = TRUE)))%>% 
  rename_at(vars(meanTemp:lGPP), function(x) paste0(x,".GMean"))

metG <- rbind(metG0, metG0, metG0, metG0)
metG$stream <- c("st11U", "st18", "st6", "st9")

Now, we log transform Q, GPP, and light then center on mean logged value (you can’t center then log).

met <- met2 %>% 
  full_join(metS, by = "stream") %>% 
  full_join(metG, by = "stream") %>% 
  #these are centered on the stream mean
  mutate(Tanom.iktCs = invKT.C - invKT.C.StMean, 
         Tanom.iktCg = invKT.C - invKT.C.GMean, 
         Tanom.Cs = meanTemp - meanTemp.StMean,
         Tanom.Cg = meanTemp - meanTemp.GMean,
         LightPerDayCs = LightPerDay - LightPerDay.StMean,
         LightPerDayCg = LightPerDay - LightPerDay.GMean,
         meanQdsCs = meanQds - meanQds.StMean,
         meanQdsCg = meanQds - meanQds.GMean,
         LightPerDay.l = log(LightPerDay),
         meanQds.l = log(meanQds),
         LightPerDay.lcs = LightPerDay.l - log(LightPerDay.StMean),
         LightPerDay.lcg = LightPerDay.l - log(LightPerDay.GMean),
         meanQds.lcs = meanQds.l - log(meanQds.StMean),
         meanQds.lcg = meanQds.l - log(meanQds.GMean),
         lGPP.cs = lGPP - lGPP.StMean,
         lGPP.cg = lGPP - lGPP.GMean,
         streamTreat = paste0(stream," ", treatment),
         stream = as.factor(stream))

Quick visualization

ER

Focus on ER response.

Access normality: log ER - not terrible

Light

Untransformed daily integrated light

log-centered daily integrated light - stream

log-centered daily integrated light - global

Is light normal? Probably close enough.

Discharge

Untransformed mean daily Q.

log-centered mean daily Q -stream

log-centered mean daily Q - global

Is it normal. Log centered mean daily Q.

Temperature as 1/kT

Untransformed temperature

Temperature anom based on mean stream temp

Temperature anom based on mean stream temp - global

Normal? - looks awesome!

GPP

Is really important and there’s super interesting things going on between ER and GPP, but I’m going to leave out.

ggplot(met, aes(y = lER, x = lGPP.cg, color = stream, shape = treatment)) +geom_point() +facet_grid(stream~.)

Normal? - S18 has particularly long tails

ggplot(met, aes(y = log(-ER_median), x = meanQds.lcg, color = Yf)) +
  geom_point() +
  facet_grid(. ~ stream)

Ugh - Q and temp are correlated These can’t be in the same model

ggplot(met, aes(y = meanQdsCg, x = Tanom.iktCg, color = stream)) +
  geom_point() +
  facet_grid(stream ~.)

met %>% 
  group_by(stream) %>% 
  summarize(currelation = cor(meanQdsCg, Tanom.iktCg))
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 4 x 2
##   stream currelation
##   <fct>        <dbl>
## 1 st11U        0.517
## 2 st18         0.714
## 3 st6          0.694
## 4 st9          0.487

ER v. temp

ggplot(met, aes(y = lER, x = Tanom.iktCg, color = Yf, shape = stream)) +
  geom_point() 

ggplot(met, aes(y = lER, x = Tanom.iktCg, color = Yf)) +
  geom_point() +
  facet_grid(stream ~.)

ggplot(met, aes(y = lER, x = invKT.C.StMean, color = Yf, shape = stream)) +
  geom_point(position = "jitter") 

ggplot(met, aes(y = lER, x = invKT.C.GMean, color = Yf, shape = stream)) +
  geom_point(position = "jitter") 

ER models

Likely maximal model

g.er.PutGLOB0 <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream) , REML = FALSE, data = met)

g.er.PutGLOB1 <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream) + 
                     s(lGPP.cg, by = stream), REML = FALSE, data = met)

g.er.PutGLOB1a <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream), REML = FALSE, data = met)


# this is perhaps the most intuitive, but it doesn't work
# g.er.PutGLOB1 <- gamm(lER ~ s(invKT.C.StMean, by = treatment) +  s(meanQds.lc, by = stream) + 
#                      s(LightPerDay.lc, by = stream) + s(Tanom.iktC, by = stream), REML = FALSE, data = met)
model.sel(g.er.PutGLOB0$lme, g.er.PutGLOB1$lme, g.er.PutGLOB1a$lme)
## Model selection table 
##                    X             family
## g.er.PutGLOB1$lme  + gaussian(identity)
## g.er.PutGLOB0$lme  + gaussian(identity)
## g.er.PutGLOB1a$lme + gaussian(identity)
##                                                                                                                                                                                                                                   random
## g.er.PutGLOB1$lme  X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g+X.3-g.3%i%g.2%i%g.1%i%g.0%i%g+X.4-g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g+X.5-g.5%i%g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g+X.6-g.6%i%g.5%i%g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g
## g.er.PutGLOB0$lme                                                                                                                                                              X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g
## g.er.PutGLOB1a$lme                                                                                                                                                             X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g
##                    df   logLik   AICc  delta weight
## g.er.PutGLOB1$lme  23 -415.529  878.7   0.00      1
## g.er.PutGLOB0$lme  15 -610.288 1251.3 372.55      0
## g.er.PutGLOB1a$lme 15 -610.288 1251.3 372.55      0
## Models ranked by AICc(x) 
## Random terms: 
## X-g = 'Xr - 1 | g'
## X.0-g.0%i%g = 'Xr.0 - 1 | g.0 %in% g'
## X.1-g.1%i%g.0%i%g = 'Xr.1 - 1 | g.1 %in% g.0 %in% g'
## X.2-g.2%i%g.1%i%g.0%i%g = 'Xr.2 - 1 | g.2 %in% g.1 %in% g.0 %in% g'
## X.3-g.3%i%g.2%i%g.1%i%g.0%i%g = 'Xr.3 - 1 | g.3 %in% g.2 %in% g.1 %in% g.0 %in% g'
## X.4-g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g = 'Xr.4 - 1 | g.4 %in% g.3 %in% g.2 %in% g.1 %in% g.0 %in% g'
## X.5-g.5%i%g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g = 'Xr.5 - 1 | g.5 %in% g.4 %in% g.3 %in% g.2 %in% g.1 %in% g.0 %in% g'
## X.6-g.6%i%g.5%i%g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g = 'Xr.6 - 1 | g.6 %in% g.5 %in% g.4 %in% g.3 %in% g.2 %in% g.1 %in% g.0 %in% g'
summary(g.er.PutGLOB1$gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## lER ~ treatment * invKT.C.StMean + s(meanQds.lcg, by = stream) + 
##     s(lGPP.cg, by = stream)
## 
## Parametric coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         0.74379    0.07204  10.325  < 2e-16 ***
## treatmentnitrogen                   1.69449    0.06688  25.337  < 2e-16 ***
## treatmentphosphorus                 1.49621    0.05325  28.095  < 2e-16 ***
## invKT.C.StMean                     -0.57992    0.14679  -3.951 8.64e-05 ***
## treatmentnitrogen:invKT.C.StMean    0.41608    0.15170   2.743  0.00626 ** 
## treatmentphosphorus:invKT.C.StMean  0.57231    0.12917   4.431 1.10e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                              edf Ref.df       F  p-value    
## s(meanQds.lcg):streamst11U 1.000  1.000   0.953    0.329    
## s(meanQds.lcg):streamst18  5.023  5.023  12.240 1.96e-11 ***
## s(meanQds.lcg):streamst6   1.000  1.000  17.223 3.75e-05 ***
## s(meanQds.lcg):streamst9   3.214  3.214  22.691 1.71e-14 ***
## s(lGPP.cg):streamst11U     1.000  1.000   0.934    0.334    
## s(lGPP.cg):streamst18      4.015  4.015  73.572  < 2e-16 ***
## s(lGPP.cg):streamst6       1.000  1.000 123.824  < 2e-16 ***
## s(lGPP.cg):streamst9       2.710  2.710  47.278  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.853   
##   Scale est. = 0.18624   n = 677
g.er.PutGLOB0.b <- getViz(g.er.PutGLOB1$gam)
check(g.er.PutGLOB0.b)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

So autocorrelation

# https://www.kaggle.com/timib1203/managing-many-models-with-tidyr-purrr-and-broom
met.PutGLOB <- met %>% 
  mutate(g.erRes = resid(g.er.PutGLOB1$lme, type = "normalized")) %>% 
  select(streamTreat, g.erRes) %>% 
  group_by(streamTreat) %>% 
  nest() %>% 
  mutate(acf = map(.x = data, .f = function(d) ggAcf(d$g.erRes) + labs(title = streamTreat)))
grid.arrange(grobs = met.PutGLOB$acf, ncol = 3)

Evaluate different corr structures

g.er.AC_1.0 <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = corARMA(p = 1,q = 0, form = ~ doy | streamTreat), REML = FALSE, data = met)

g.er.AC_1.1 <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = corARMA(p = 1,q = 1, form = ~ doy | streamTreat), REML = FALSE, data = met)

g.er.AC_1.2 <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = corARMA(p = 1,q = 2, form = ~ doy | streamTreat), REML = FALSE, data = met)

g.er.AC_2.0 <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = corARMA(p = 2,q = 0, form = ~ doy | streamTreat), REML = FALSE, data = met)


g.er.AC_2.1 <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = corARMA(p = 2,q = 1, form = ~ doy | streamTreat), REML = FALSE, data = met)

g.er.AC_2.2 <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = corARMA(p = 2,q = 2, form = ~ doy | streamTreat), REML = FALSE, data = met)
# doesn't work
g.er.AC_3.0 <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = corARMA(p = 3,q = 0, form = ~ doy | streamTreat), REML = FALSE, data = met)

g.er.AC_3.1 <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = corARMA(p = 3,q = 1, form = ~ doy | streamTreat), REML = FALSE, data = met)
# 
g.er.AC_3.2 <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = corARMA(p = 3,q = 2, form = ~ doy | streamTreat), REML = FALSE, data = met)

Model selection

model.sel(g.er.AC_1.0$lme, g.er.AC_1.1$lme, g.er.AC_1.2$lme, g.er.AC_2.0$lme,
          g.er.AC_2.1$lme, g.er.AC_2.2$lme, g.er.AC_3.0$lme, g.er.AC_3.1$lme, g.er.AC_3.2$lme) 
## Model selection table 
##                 X             family df  logLik  AICc delta weight
## g.er.AC_3.2$lme + gaussian(identity) 20 -35.365 112.0  0.00  0.933
## g.er.AC_3.0$lme + gaussian(identity) 18 -41.171 119.4  7.37  0.023
## g.er.AC_3.1$lme + gaussian(identity) 19 -40.472 120.1  8.09  0.016
## g.er.AC_1.2$lme + gaussian(identity) 18 -41.761 120.6  8.55  0.013
## g.er.AC_2.2$lme + gaussian(identity) 19 -40.709 120.6  8.56  0.013
## g.er.AC_1.1$lme + gaussian(identity) 17 -45.297 125.5 13.51  0.001
## g.er.AC_2.0$lme + gaussian(identity) 17 -46.246 127.4 15.41  0.000
## g.er.AC_1.0$lme + gaussian(identity) 16 -48.863 130.5 18.54  0.000
## g.er.AC_2.1$lme + gaussian(identity) 18 -48.438 133.9 21.91  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
## 'Xr - 1 | g', 'Xr.0 - 1 | g.0 %in% g', 'Xr.1 - 1 | g.1 %in% g.0 %in% g', 'Xr.2 - 1 | g.2 %in% g.1 %in% g.0 %in% g'

Much better. S6p, s6n, s9a, s9p, s11u still have some minor problems

# https://www.kaggle.com/timib1203/managing-many-models-with-tidyr-purrr-and-broom
# https://jroy042.github.io/nonlinear/week4.html
# https://stackoverflow.com/questions/47586110/autocorrelation-in-generalized-additive-models-gam
g.er.AC_best.acf <- met %>% 
  mutate(g.erRes = resid(g.er.AC_3.2$lme, type = "normalized")) %>% 
  select(streamTreat, g.erRes) %>% 
  group_by(streamTreat) %>% 
  nest() %>% 
  mutate(acf = map(.x = data, .f = function(d) ggAcf(d$g.erRes) + labs(title = streamTreat)))
grid.arrange(grobs = g.er.AC_best.acf$acf, ncol = 3)

Let’s do some model selection. Thinking about this as our base model has Q, temp, and light. I’m not going to obsess about selecting those. Instead, I’ll just focus on Annual Temp and treatment

going to use a simpler one, problems fitting models with 3.0

BestERCor <- corARMA(p = 3,q = 2, form = ~ doy | streamTreat)
# did quick check to see if I need T anom by stream - yes massive improvement in AIC
# also did check to see if I needed GPP x stream - nope! 
# GLOBAL CENTERED Q
g.er.NutXT <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.NutpT <- gamm(lER ~ treatment + invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)
g.er.NutpTs <- gamm(lER ~ treatment + invKT.C.StMean +  s(meanQds.lcs, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.NutpTb <- gamm(lER ~ treatment + invKT.C.StMean +  s(meanQds.lcg),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.NutXst <- gamm(lER ~ treatment*stream +  s(meanQds.lcg, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.Nutpst <- gamm(lER ~ treatment + stream +  s(meanQds.lcg, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.Nut <- gamm(lER ~ treatment +  s(meanQds.lcg, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.T <- gamm(lER ~ invKT.C.StMean +  s(meanQds.lcg, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.st <- gamm(lER ~ stream +  s(meanQds.lcg, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.null <- gamm(lER ~ s(meanQds.lcg, by = stream),
                correlation = BestERCor, REML = FALSE, data = met) 

# STREAM CENTERED Q
g.er.NutXTs <- gamm(lER ~ treatment*invKT.C.StMean +  s(meanQds.lcs, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.NutpTs <- gamm(lER ~ treatment + invKT.C.StMean +  s(meanQds.lcs, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.NutXsts <- gamm(lER ~ treatment*stream +  s(meanQds.lcs, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.Nutpsts <- gamm(lER ~ treatment + stream +  s(meanQds.lcs, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.Nuts <- gamm(lER ~ treatment +  s(meanQds.lcs, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.Ts <- gamm(lER ~ invKT.C.StMean +  s(meanQds.lcs, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.sts <- gamm(lER ~ stream +  s(meanQds.lcs, by = stream),
                correlation = BestERCor, REML = FALSE, data = met)

g.er.nulls <- gamm(lER ~ s(meanQds.lcs, by = stream),
                correlation = BestERCor, REML = FALSE, data = met) 

Probably a tie between treatment and treatment + temp But the model with a stream centered Q is wicked better

model.sel(g.er.NutXT$lme, g.er.NutpT$lme, g.er.NutXst$lme, g.er.Nutpst$lme, g.er.Nut$lme, g.er.T$lme, g.er.st$lme, g.er.null$lme, g.er.NutpTb$lme,
          g.er.NutXTs$lme, g.er.NutpTs$lme, g.er.NutXsts$lme, g.er.Nutpsts$lme, g.er.Nuts$lme, g.er.Ts$lme, g.er.sts$lme, g.er.nulls$lme)
## Model selection table 
##                  X             family
## g.er.NutpTs$lme  + gaussian(identity)
## g.er.Nutpsts$lme + gaussian(identity)
## g.er.NutXsts$lme + gaussian(identity)
## g.er.NutXTs$lme  + gaussian(identity)
## g.er.Nuts$lme    + gaussian(identity)
## g.er.Nut$lme     + gaussian(identity)
## g.er.NutpT$lme   + gaussian(identity)
## g.er.Nutpst$lme  + gaussian(identity)
## g.er.NutXT$lme   + gaussian(identity)
## g.er.NutXst$lme  + gaussian(identity)
## g.er.nulls$lme   + gaussian(identity)
## g.er.sts$lme     + gaussian(identity)
## g.er.null$lme    + gaussian(identity)
## g.er.T$lme       + gaussian(identity)
## g.er.Ts$lme      + gaussian(identity)
## g.er.st$lme      + gaussian(identity)
## g.er.NutpTb$lme  + gaussian(identity)
##                                                                     random df
## g.er.NutpTs$lme  X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 18
## g.er.Nutpsts$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 20
## g.er.NutXsts$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 26
## g.er.NutXTs$lme  X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 20
## g.er.Nuts$lme    X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 17
## g.er.Nut$lme     X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 17
## g.er.NutpT$lme   X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 18
## g.er.Nutpst$lme  X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 20
## g.er.NutXT$lme   X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 20
## g.er.NutXst$lme  X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 26
## g.er.nulls$lme   X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 15
## g.er.sts$lme     X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 18
## g.er.null$lme    X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 15
## g.er.T$lme       X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 16
## g.er.Ts$lme      X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 16
## g.er.st$lme      X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 18
## g.er.NutpTb$lme                                                        X-g 12
##                   logLik  AICc delta weight
## g.er.NutpTs$lme  -32.672 102.4  0.00  0.447
## g.er.Nutpsts$lme -31.127 103.5  1.15  0.251
## g.er.NutXsts$lme -25.728 105.6  3.23  0.089
## g.er.NutXTs$lme  -32.594 106.5  4.08  0.058
## g.er.Nuts$lme    -35.776 106.5  4.10  0.058
## g.er.Nut$lme     -35.840 106.6  4.22  0.054
## g.er.NutpT$lme   -35.463 108.0  5.58  0.027
## g.er.Nutpst$lme  -34.297 109.9  7.49  0.011
## g.er.NutXT$lme   -35.365 112.0  9.63  0.004
## g.er.NutXst$lme  -29.408 113.0 10.59  0.002
## g.er.nulls$lme   -43.347 117.4 15.04  0.000
## g.er.sts$lme     -42.082 121.2 18.82  0.000
## g.er.null$lme    -45.494 121.7 19.33  0.000
## g.er.T$lme       -45.321 123.5 21.08  0.000
## g.er.Ts$lme      -46.578 126.0 23.60  0.000
## g.er.st$lme      -44.544 126.1 23.74  0.000
## g.er.NutpTb$lme  -59.847 144.2 41.78  0.000
## Models ranked by AICc(x) 
## Random terms: 
## X-g = 'Xr - 1 | g'
## X.0-g.0%i%g = 'Xr.0 - 1 | g.0 %in% g'
## X.1-g.1%i%g.0%i%g = 'Xr.1 - 1 | g.1 %in% g.0 %in% g'
## X.2-g.2%i%g.1%i%g.0%i%g = 'Xr.2 - 1 | g.2 %in% g.1 %in% g.0 %in% g'
summary(g.er.NutpTs$gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## lER ~ treatment + invKT.C.StMean + s(meanQds.lcs, by = stream)
## 
## Parametric coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           0.8139     0.2218   3.669 0.000263 ***
## treatmentnitrogen     1.7862     0.3011   5.933 4.81e-09 ***
## treatmentphosphorus   1.3471     0.2994   4.500 8.05e-06 ***
## invKT.C.StMean       -1.0141     0.3709  -2.734 0.006426 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                              edf Ref.df      F  p-value    
## s(meanQds.lcs):streamst11U 2.484  2.484  7.461 0.000658 ***
## s(meanQds.lcs):streamst18  3.597  3.597  4.288 0.002139 ** 
## s(meanQds.lcs):streamst6   4.039  4.039 11.987 1.71e-09 ***
## s(meanQds.lcs):streamst9   5.030  5.030 20.046  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.609   
##   Scale est. = 0.46633   n = 677
g.er.NutpT.b <- getViz(g.er.NutpTs$gam)

This model violates some shit!!!

check(g.er.NutpT.b)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

hist(resid(g.er.NutpT$lme, type = "normalized"))

g.er.NutpT.acf <- met %>% 
  mutate(g.erRes = resid(g.er.NutpTs$lme, type = "normalized")) %>% 
  select(streamTreat, g.erRes) %>% 
  group_by(streamTreat) %>% 
  nest() %>% 
  mutate(acf = map(.x = data, .f = function(d) ggAcf(d$g.erRes) + labs(title = streamTreat)))
grid.arrange(grobs = g.er.NutpT.acf$acf, ncol = 3)

met$BestMod <- predict.gam(g.er.NutpTs$gam, met, type = "response")
ggplot(met, aes(y = lER, x = BestMod, color = treatment)) +
  geom_point()+
  facet_grid(stream ~.) +
  geom_abline(intercept = 0, slope =1)

print(plot(g.er.NutpT.b, allTerms = T), pages = 1)

In the splines, this shows all the data!

blah <- plot(g.er.NutpT.b, allTerms = 1) +
  l_points(color = "blue") + l_fitLine(linetype = 3) + l_fitContour() + 
      l_ciLine(colour = 2) + l_ciBar() + l_fitPoints(size = 1, col = 2) + theme_get() + labs(title = NULL)

print(blah, pages = 1)

save.image("01_ERanalysis_rdat")
# load("01_ERanalysis_rdat")